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G\ . ABSTRACT 



We present numerical simulations of off-center collisions between galaxy clusters made 
using a new hydrodynamical code based on the piecewise-parabolic method (PPM) and 
an isolated multigrid potential solver. The current simulations follow only the intraclus- 
ter gas. We have performed three high-resolution (256 x 128 2 ) simulations of collisions 
between equal-mass clusters using a nonuniform grid with different values of the impact 
parameter (0, 5, and 10 times the cluster core radius). Using these simulations we have 
studied the variation in equilibration time, luminosity enhancement during the collision, 
jy-) | and structure of the merger remnant with varying impact parameter. 

We find that in off-center collisions the cluster cores (the inner regions where the 
pressure exceeds the ram pressure) behave quite differently from the clusters' outer 
regions. A strong, roughly ellipsoidal shock front, similar to that noted in previous 
simulations of head-on collisions, enables the cores to become bound to each other by 
dissipating their kinetic energy as heat in the surrounding gas. These cores survive well 
into the collision, dissipating their orbital angular momentum via spiral bow shocks. 
After the ellipsoidal shock has passed well outside the interaction region, the material 
left in its wake falls back onto the merger remnant formed through the inspiral of the 
cluster cores, creating a roughly spherical accretion shock. For less than one-half of a 
sound crossing time after the cores first interact the total X-ray luminosity increases 
by a large factor; the magnitude of this increase depends sensitively on the size of the 
impact parameter. 

Observational evidence of the ongoing collision, in the form of bimodality and dis- 
tortion in projected X-ray surface brightness and temperature maps, is present for 1-2 
sound crossing times after the collision, but only for special viewing angles. The remnant 
actually requires at least five crossing times to reach virial equilibrium. Since the sound 
crossing time can be as large as 1-2 Gyr, the equilibration time can thus be a substantial 
fraction of the age of the universe. The final merger remnant is very similar for impact 
parameters of zero and five core radii. It possesses a roughly isothermal core, with cen- 
tral density and temperature twice the initial values for the colliding clusters. Outside 
the core the temperature drops as r _1 , and the density roughly as r -3 8 . The core radius 
shows a small increase due to shock heating during the merger. For an impact parameter 
of ten core radii the core of the remnant possesses a more flattened density profile, with 
a steeper dropoff outside the core. In both off-center cases the merger remnant rotates, 
but only for the ten-core-radius case does this appear to have an effect on the structure 
of the remnant. 
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1. Introduction 

Galaxy clusters are the largest bound structures in 
the present-day universe. However, they are not the 
oldest. For galaxy clusters the notion of 'age' does not 
necessarily refer to the time elapsed since a specific 
formation event occurred, after which they remained 
either static in their properties or very slowly evolv- 
ing. Indeed, a growing body of observational and the- 
oretical evidence indicates that for many clusters 'for- 
mation' means an ongoing sequence of mergers and 
interactions with other clusters. Understanding this 
process is important not only from the standpoint of 
understanding cluster evolution, but also because of 
its implications for cosmology. 

Recent ROSAT and ASCA observations of rich 
clusters of galaxies have shown that many of them are 
either in the process of merging with another clus- 
ter or have recently undergone such a merger. The 
ROSAT results, because of poor spectral resolution 
but good spatial resolution, emphasize the variation 
in gas density and hence the depth of the potential 
well in a cluster. These include observations of the 
Coma cluster (White, Briel, & Henry 1993), Abell 
2256 (Briel et al. 1991), and Abell 2163 (Elbaz, Ar- 
naud, & Bohringer 1995). The signatures of merging 
events in these observations include distorted X-ray 
isophotcs (presumably elongated along the collision 
axes); offsets of the gas centroid, as indicated by the 
X-ray emission, and the center of mass of the galaxy 
distribution; and bimodal distributions of emission. 
The ROSAT results are consistent with earlier results 
from Einstein observations which suggested that at 
least 22% of clusters show significant X-ray substruc- 
ture (Forman & Jones 1994). In the more recent 
ASCA observations of Abell 754 (Henriksen & Marke- 
vitch 1996), which emphasize the variation in gas tem- 
perature through its influence on the X-ray spectrum, 
the presence of a ridge of high-temperature emission 
has been taken as indicative of an ongoing merger 
event. However, significant temperature variations 
have not been seen by ASCA in other merger can- 
didates identified from ROSAT data, including Abell 
2256 and Abell 2163 (Markevitch 1996). 

Although we do not yet completely understand the 
physical mechanisms behind and cosmological impli- 
cations of galaxy cluster mergers, it is not surpris- 
ing that they should occur frequently. In hierarchical 
models of large-scale structure formation, objects like 
galaxy clusters form through the merging and then ac- 



cretion of smaller clumps of matter. In these models 
the present-day fraction of clusters which have sig- 
nificant substructure can give us constraints on the 
value of the density parameter fl (Richstone, Loeb, & 
Turner 1992), and the mass density profiles resulting 
from mergers provide information about the spectrum 
of primordial density fluctuations (Hoffman & Sha- 
ham 1985). However, although analytical approaches 
based on the spherical top-hat model and the Press- 
Schechter (1974) mass function have provided much 
insight into these implications, the complex physical 
interactions present during a cluster merger can only 
be studied adequately through numerical simulation. 

To date simulations of galaxy cluster evolution (ei- 
ther A-body or mesh-based) have typically followed 
one of three approaches. Some simulations have fo- 
cused on the collapse and relaxation of a single over- 
dense region which has been prepared as a 'typical' 
example given a large-scale structure model (e. g. 
Evrard 1990; Evrard, Metzler, & Navarro 1996). Oth- 
ers have produced clusters within large-scale structure 
simulations, then followed their evolution separately 
from the rest with greater spatial resolution (e. g. 
Katz & White 1993). Both types of simulation pro- 
duce realistic clusters in an environment which per- 
mits comparison of large-scale structure theories but 
which renders analysis of individual merging events 
difficult because of the uncontrolled initial conditions 
for each event and the frequent occurrence of multiple 
simultaneous mergers. 

A third type of cluster evolution simulation has 
focused on collisions of idealized clusters or subclus- 
ters (e. g. Roettiger, Burns, & Loken 1993; Schindlcr 
& Midler 1993; Pearcc, Thomas, & Couchman 1994; 
Roettiger, Stone, & Mushotzky 1997) to gain in- 
sight into the merger process without complicating 
the problem with multiple mergers. This approach 
also permits mergers to be studied with better reso- 
lution than is possible in large-scale structure simula- 
tions. 

Roettiger et al. (1993) simulated the head-on colli- 
sion of two King- model clusters of unequal mass using 
the ZEUS-3D finite-difference code with a nonuniform 
grid for the gas dynamics and the Hernquist (1987) 
treecode for the dark matter. The resolution of these 
simulations was such that the smallest zones were one- 
half the core radius of the smaller cluster. During the 
collision they observed a double-peaked X-ray emis- 
sion profile and a bar of X-ray emission perpendicular 
to the collision axis which resulted from a strong cen- 
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tral shock (Mach number ~ 4). They also observed 
a high-velocity (~ 2800 km/s), ordered flow of gas 
through the center of the dominant cluster which they 
suggested might disrupt cooling flows and be respon- 
sible for the bending of radio jets in cluster galaxies. 
Although their simulation did not include cooling, a 
later simulation by Burns et al. (1994) of a collision 
between the Coma cluster and the NGC 4839 group, 
using the same code but including cooling and hav- 
ing somewhat higher core resolution, obtained similar 
results. 

Schindler and Miiller (1993) also studied head- 
on collisions derived from slightly less controlled ini- 
tial conditions using the piecewise-parabolic method 
(PPM) for the gas dynamics and Aarseth's (1972) 
code for collisionless galaxy halos. During the course 
of their collisions they observed multiple interacting 
shocks which led to X-ray luminosity enhancements of 
about a factor of 1.5 and temperature enhancements 
of a factor of 5, with a non- isothermal final state as 
the result. These calculations used a uniform grid 
with about 2.5 cells per core radius. 

Neither Rocttigcr ct al. (1993), Burns et al., nor 
Schindler and Miiller included the baryonic contri- 
bution to the gravitational potential, thus limiting 
their conclusions to gas-poor clusters (less than 5%- 
10% of the total mass). However, some clusters can 
have as much as 30% (for a Hubble constant of 50 
km s" 1 Mpc -1 ) of their total mass in the intraclus- 
ter medium (David, Jones, & Forman 1995; White 
& Fabian 1995). The virial estimates of total clus- 
ter masses on which this ROSAT finding is based are 
supported by recent simulation work by Evrard et al. 
(1996). The gas in these clusters may therefore have 
a significant effect on the potential, rendering a self- 
consistent treatment necessary. 

Pearce et al. (1994) performed head-on collision 
simulations of equal-mass clusters using an adaptive 
smoothed-particlc hydrodynamics (SPH) code which 
does include both gas and dark matter in the po- 
tential. They found that the gas, which made up 
1/8 of their clusters' total mass, tends to produce 
an extended constant-density core after the collision, 
whereas the dark matter component does not. How- 
ever, while they note the presence of shocks due to 
the collision of the gas components of the clusters, 
they do not discuss the observational appearance of 
the merging system in detail or compare it to the 
finite-difference results. More recently, Roettiger et 
al. (1997) have used a PPM/particle- mesh code to 



simulate a collision intended to reproduce the ob- 
servations of A754, varying collision parameters at 
low resolution and then resimulating with the most 
promising parameter values at high resolution. 

With the exception of Roettiger et al. (1997), none 
of these authors have considered off-center collisions 
in detail, although Pearce et al. mention that pre- 
liminary results from such simulations indicate that 
the gas receives additional support from rotation in 
the merger remnant. Interest in off-center collisions 
is growing, however, as ASCA temperature maps of 
clusters such as A754 (Henriksen & Markevitch 1996) 
begin to give us more detailed information about off- 
center cluster mergers in progress. The net angular 
momentum imparted to a cluster by many such col- 
lisions is probably small, since the angular momenta 
brought by collisions in different directions will tend 
to cancel each other out. Also, the net spins in- 
duced in collapsed objects by tidal torques due to 
inhomogencities in the large-scale mass distribution 
have been shown to be small (Peebles 1993). However, 
during their pre-collision evolution, colliding subclus- 
ters should experience some deflection due to distant 
collisions with other subclusters. Given the great dif- 
ference in behavior to be expected between the limit- 
ing cases of a simple distant deflection and a head-on 
collision, we would do well to consider the dynamics 
of intermediate cases and its effects on the long-term 
evolution and observable properties of clusters. 

In order to study galaxy cluster evolution, and in 
particular off-center collisions of galaxy clusters, we 
have developed a parallel hydrodynamics code based 
on PPM which self-consistently solves for the gravi- 
tational potential of the gas using a multigrid algo- 
rithm with isolated boundary conditions. The code 
uses a nonuniform grid, permitting it to resolve the 
inner regions of the colliding clusters while simultane- 
ously tracking the development of the shocks which 
result from the collision. We have added radiative 
cooling and linked the PPM code to a particle-mesh 
iV-body solver for dark matter (written by Scott Do- 
delson), but the simulations described in this paper 
include only nonradiating gas. While these gas-only 
calculations necessarily exclude the substantial con- 
tribution of collisionless effects to cluster evolution, 
they nevertheless represent an important benchmark 
for more realistic calculations including dark matter, 
because it is far more straightforward to understand 
both the physical and numerical properties of the hy- 
drodynamical results in isolation. This is analogous 
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to the approach taken prior to the advent of com- 
bined TV-body/hydro codes, in which the collisionless 
dynamics of the dark matter was studied in isolation 
using iV-body codes with well-understood properties. 
In a forthcoming paper we will discuss simulations 
performed with our new hybrid code, using the re- 
sults presented here as a basis of comparison. 

The paper is organized as follows. Section 2 gives 
a brief description of the code and the initial and 
boundary conditions used; a more complete descrip- 
tion, including the equations used and the results of 
running the code on several standard test problems, 
will be presented in another forthcoming paper. Sec- 
tion 3 describes a resolution study performed using a 
single cluster at rest at the center of the grid. Section 
4 describes results from three high-resolution simu- 
lations of collisions between equal-mass clusters at 
different impact parameters. Section 5 discusses the 
equilibration times and X-ray brightening observed 
in these calculations, as well as the thermal and ro- 
tational properties of the merger remnants. Section 6 
summarizes our conclusions. 

2. Numerical methods 
2.1. Hydrodynamics 

To solve the equations of hydrodynamics for the 
cluster gas we use the direct Eulerian formulation 
of the piecewise-parabolic method (PPM; Colella & 
Woodward 1984) on a nonuniform Cartesian grid. In 
contrast to most numerical studies of cluster evolution 
which use smoothed particle hydrodynamics (SPH) 
(e. g. Evrard 1990) for its ability to dynamically ad- 
just spatial resolution, simulations using PPM offer 
better shock handling for comparable spatial resolu- 
tion, since they require very little artificial viscos- 
ity. SPH resolves shocks only poorly and requires 
substantial amounts of artificial viscosity. Further- 
more, in low-density regions where SPH represents 
the gas with relatively few particles, PPM obtains 
more accurate temperatures (Kang et al. 1994). The 
error properties of finite- volume methods in general, 
of which PPM is a specific example, are better under- 
stood than those of SPH. For these reasons we believe 
PPM to be the more appropriate method at present 
for studying the evolution of the collisional compo- 
nent of galaxy clusters. 

Others have used PPM to study cluster evolution. 
Schindler and Miiller (1993) used PPM together with 
the iV-body solver of Aarseth (1972) to evolve the gas 



in their simulation on a 60 3 uniform mesh covering 
little more than the two colliding clusters they stud- 
ied. Bryan et al. (1994) used PPM on a much larger 
mesh with 270 3 cells to study the evolution of many 
poorly resolved clusters over a volume large enough to 
contain a statistical sample of clusters. Roettiger et 
al. (1993) used a non-Godunov-based finite-difference 
method similar to PPM in its use of parabolic in- 
terpolants to follow two merging clusters; their grid 
contained 187 x 65 2 cells and had nonuniform spac- 
ing which permitted good resolution of the cluster 
cores. More recently, Roettiger et al. (1997) have 
used a PPM/particle-mesh hybrid code to simulate 
the evolution of Abell 754 on a 512 x 256 2 grid with a 
maximum resolution of 1/10 the cluster core radius. 

However, with the exception of Bryan et al. and 
Roettiger et al. (1997), none of these authors self- 
consistently solved for the contribution of the gas to 
the potential. This is a serious shortcoming, one not 
shared by the SPH codes; even a gas mass fraction 
of 10% can significantly affect the potential on small 
scales if the gas and dark matter distributions are 
not the same. In our code we therefore include the 
gas density in the computation of the potential, which 
is described in the next subsection. Although the re- 
sults presented here do not include dark matter, we 
have linked our hydrodynamical code to an iV-body 
code based on the particle-mesh method, with the two 
codes sharing the same potential solver. We plan to 
report on results obtained with the new hybrid code 
in the near future. The new code also includes radia- 
tive cooling, which we do not include in the results 
presented here. 

In order to look for observational signatures of re- 
cent on- or off-center merger events, one must be 
able to study a computational volume which is large 
enough to contain the shocks produced by the merger 
until they cease to be important. Furthermore, our 
single-cluster test runs, described in the next section, 
suggest that to do the problem correctly one must use 
a grid spacing in the inner regions of the clusters of at 
most one-half the core radius r c . In order to begin to 
reconcile these two requirements we have used for our 
high- resolution runs a nonuniform grid of 256 x 128 2 
zones. The innermost A x i nnor = 192 zones of the 
ai-axis have uniform widths equal to 0.28r c . The in- 
nermost -/V yZj i nncr = 72 zones in y and z have uniform 
widths equal to 0.39r c . Outside the uniformly gridded 
region in each coordinate direction the zone widths in- 
crease at the rate of 5% per zone as one moves away 
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from the center of the grid. The size of the uniformly 
gridded region was chosen to maximize resolution and 
minimize any effects due to grid nonuniformity in the 
region where most of the complex shocks form during 
the simulation. In terms of the initial cluster radius 
R = 9r c , our computational grid has physical dimen- 
sions 11R x 8.1R x S.1R. 

In test calculations of the Sedov- Taylor explosion 
problem (Sedov 1959) we have found that using first- 
order operator splitting to create a three-dimensional 
PPM scheme does not adequately maintain the 90- 
degree rotational symmetry of a uniform grid, partic- 
ularly in the central part of the grid where a rar- 
efied region develops in the Sedov- Taylor problem. 
We therefore use a second-order splitting, generalized 
from Strang (1968) to three dimensions. Although 
this doubles the cost of a hydrodynamical timestep 
relative to first-order splitting, the spherically sym- 
metric, centrally concentrated nature of our initial 
density and pressure field makes it important to pre- 
serve this symmetry. 

In a realistic simulation of cluster evolution, tides 
from mass concentrations outside of the grid would 
affect the simulated objects. Mass inflow through the 
grid boundaries would also be expected. However, 
in our simulations we are studying the behavior of a 
single merging system, so vacuum boundaries are ap- 
propriate. To implement these we use the standard 
zero-gradient boundaries with a slight difference. If 
purely zero-gradient boundaries are used, any inward- 
directed velocity in the interior zones adjacent to the 
boundary will result in a completely artificial inward 
mass flux. Since this mass flux adds to the depth 
of the potential, adding to the inward acceleration in 
these zones, this problem can destabilize the entire 
calculation. To prevent such artificial inward flows 
we modify the zero-gradient prescription when veloc- 
ities in the interior zones adjacent to the boundaries 
are directed inward. For these zones we apply Dirich- 
let boundary conditions to the component of veloc- 
ity normal to the boundary. All other variables, in- 
cluding velocity components parallel to the boundary, 
are given zero-gradient boundary conditions as usual. 
This suffices for the interpolation step in PPM; for the 
flux computation step we set the boundary flux of all 
conserved quantities (mass, momentum, and energy) 
to zero if it would otherwise be directed inward. 



2.2. Gravitation 

At each timestep of the hydrodynamical code we 
solve the Poisson equation for the gravitational po- 
tential using a full multigrid/nested iteration code 
(Brandt 1977; Hackbusch 1985). This method ac- 
celerates the convergence of standard elliptic relax- 
ation techniques by iteratively solving on grids of in- 
creasing coarseness. Because long-wavelength errors 
are responsible for the slow convergence of standard 
iteration techniques, solving on coarse grids acceler- 
ates convergence by turning long-wavelength modes 
into short-wavelength modes. The resulting speedup 
is sufficient to make multigrid methods competitive 
with direct methods, such as those based on the fast 
Fourier transform, even for nonlinear equations or un- 
usual boundary conditions. 

In our case the primary reasons for using a multi- 
grid method are the need for ease and portability of 
parallelization and the nature of the boundary condi- 
tions for our problem. FFT-based Poisson solvers are 
most efficient when working with periodic, Dirichlet, 
or Neumann boundaries on a uniform grid. However, 
in the colliding cluster problem, it is more natural to 
permit the potential to go to zero at infinity. Meth- 
ods exist to handle these 'isolated' boundaries using 
direct solvers (e. g. James 1977) but suffer from the 
disadvantages of the direct methods on which they are 
based. To handle isolated boundaries with multigrid 
we have developed a method similar to that of James. 
We first obtain a solution with Dirichlet boundaries, 
then use this to compute image charges on the bound- 
aries under the assumption that the potential is zero 
not only on the boundary but everywhere outside it 
as well. We then compute the isolated potential of 
the image charges and subtract it from the Dirich- 
let solution. To calculate the isolated potential of the 
image charges, we use a truncated spherical harmonic 
expansion of the image charge distribution about its 
center of mass to estimate boundary values for the po- 
tential, then solve Laplace's equation in the interior 
using these boundary values. Although this method 
requires two calls to the multigrid solver, we only need 
to compute image moments and boundary values for 
points on the boundary. We have tested this potential 
solver with mass distributions drawn from collision 
runs performed on a nonuniform 128 x 64 2 grid simi- 
lar to that described in the previous section, and we 
find that an acceptable compromise between compu- 
tation time and accuracy is reached if one truncates 
the multipole expansion at I = 10. The typical change 
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in the potential <j> in going to larger i is less than 1%. 

2.3. Initial conditions 

We are interested in gaining insight into the physics 
of cluster collisions without the complication of initial 
conditions drawn from a large-scale structure simula- 
tion. By examining the parametric dependence of the 
observed features of such events we can later check 
cosmological theories by determining what distribu- 
tion of collision parameters they produce. Accord- 
ingly we have chosen for our initial conditions nonsin- 
gular isothermal spheres, approximated by the modi- 
fied Hubble law, 



p(r) 



1 + (r/r c y 



3/2 ■ 



(1) 



where po is the central mass density and r c is the 
cluster core radius, as in the paper by Pearce et al. 
(1994). However, instead of truncating the distribu- 
tion (1) beyond the radius R = 16r c as they did, we 
have used the smaller value of 9r c . The exact value 
of R/r c , if it is of order 10 or so, is not important for 
the dynamics of the inner regions of the clusters, but 
it does affect our ability to resolve the cluster cores 
properly. We assume that the density perturbations 
which created the clusters have entered the strongly 
nonlinear regime well before the start of the simu- 
lation, so we use a background spacetime which is 
flat and nonexpanding. Since PPM is a finite-volume 
scheme which uses cell-averaged quantities (in con- 
trast to a finite-difference scheme, which manipulates 
the values of fields at grid points), we set the density 
in each cell using Monte Carlo averaging of the den- 
sity profile (1). The temperature is specified by the 
core radius and central number density via the stan- 
dard definition of the core radius (Binney & Tremaine 
1987): 

fc^init 4 2 / \ 

— = r G Po r c . (2) 

We set the pressure in each cell by Monte Carlo- 
averaging p(r)T together with the density. To convert 
temperature to pressure we use 



P = 



pkT-u 



p 



(3) 



where k is Boltzmann's constant. Here p is different 
from the average atomic mass (m) because of the con- 
tribution of electrons to the gas pressure; assuming a 



primordial composition (25% helium, 75% hydrogen 
by mass), we have p « 0.59mn and (m) w 1.2toh. 
As our gas is hot enough to be completely ionized, we 
use 7 = 5/3 as our ratio of specific heats. 

In the simulations presented here, both clusters 
have the same central density and core radius, hence 
the same mass and temperature. We use units in 
which each cluster initially has unit mass, unit ra- 
dius, and unit temperature. The time unit is one 
sound crossing time, 



= — = 9r c 



p 



7 fcT ir 



(4) 



so in our units the initial sound speed c s is also unity. 
In these units the central density po is approximately 
30.5, and the core radius r c is 1/9. The scaling of 
our simulations to physical units is determined by our 
choices for the values of mn/k and Newton's gravita- 
tional constant G; we use mn/k sa 2.86i 2 c Ti n i t i? -2 
and G w 1.14i? 3 M _1 i~ 2 . Outside the two clusters, 
the density is set to 10~ 8 p(R), and the temperature 
is the same as it is inside the clusters. 

The two clusters are initialized with centers in the 
rry-plane, separated by a distance 6i n ;t (the impact 
parameter) in the y-direction and a distance d\mt in 
the x-direction (Figure 1). For the calculations de- 
scribed here we performed three 256 x 128 x 128 runs 
using different values of the impact parameter &i n ;t. 
In run A, &i„i t = 0; in run B, 6; n i t = 5r c ; and in run 
C, 6; n ;t = 10r c . The value of rfi n it was held fixed at 4i? 
for all of the runs, so that the clusters started close 
enough together to collide within a short time but far 
enough apart so that any numerical artifacts due to 
the initial conditions would be smoothed out before 
the cores interacted. The initial relative velocity, v- m \.t , 




Fig. 1. — Initial conditions for cluster collision runs. 
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was chosen to be approximately equal to the free-fall 
velocity for the two clusters at a distance rf; n i t , and 
was directed along the x-axis. All three runs used the 
same value for Vi n a (unity) . The results of these runs 
are discussed in Sections 4 and 5. 



3. Single-cluster tests 



To study the behavior of our code on a simple yet 
realistic test problem, wc simulated a single cluster at 
the center of a nonuniform grid similar to that used in 
our collision simulations. Since the gas density profile 
is approximately that of an isothermal sphere, these 
initial conditions should produce an approximately 
static configuration. We performed three runs at grid 
sizes of 32 3 (run Tl), 64 3 (run T2), and 128 3 (run 
T3) using nonuniform gridding in all three dimen- 
sions. The physical box size was kept constant at 
L « 6i?, and the grid nonuniformity, described in 
the previous section, was kept constant at 5%. This 
leaves the number of zones in the uniform inner re- 
gion, iVinnor, and the zone size in this region, A, as free 
parameters. The values of these parameters used in 
the test runs are summarized in Table 1. In each run 
the central density p and core radius r c were given 
the same values used in the collision runs. Varying 
the resolution in this way permitted us to determine 
how core resolution affects the accuracy of the solu- 
tion in the central regions of the cluster, and also to 
control for the deviation of our initial conditions from 
true isothermal spheres in the collision runs. 

Figure 2 presents the average kinetic (T), internal 
(U), and potential (W) energies of the cluster in each 
of the three test runs. For each run the total energy 
{£ = T + U + W) , less the amount lost from the grid 
between t — t sc and t — 3t sc , is constant to better 
than 1%, with the energy differing slightly between 
the runs because of slight differences in the size of 
the box used. However, in each run during the first 
2-3 sound crossing times the internal and potential 
energies decrease in magnitude, and the kinetic en- 
ergy briefly increases from zero. This occurs because 
our initial conditions contain an artificial shock where 
the Hubble-law profile is truncated. To keep the near 
vacuum outside the clusters in pressure equilibrium 
with the cluster, we would have to set the tempera- 
ture in these regions to a very large value, forcing the 
timestep to be unacceptably small (via the Courant 
condition). We therefore allow a pressure gradient to 
exist at the edges of the cluster. Initially this gradi- 
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Fig. 2. — The upper plot shows the kinetic, internal, 
and potential energies in the single cluster tests. The 
total energy is also plotted. The lower plot shows 
the time derivative of the total energy. Dashed lines 
indicate run Tl (32 3 ), dot-dashed lines, run T2 (64 3 ), 
and solid lines, run T3 (128 3 ). 
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Table 1: Grid parameters used in the single-cluster test runs. 



Run 


Grid size 


1 y inner 


A/r c 


Tl 


32 a 


16 


1.36 


T2 


64 3 


32 


0.582 


T3 


128 3 


64 


0.212 



ent drives a rapid expansion of the 'vacuum material' 
outside the cluster. The cluster itself also expands, 
cooling by about a factor of two and driving about 
15% of the total mass off the grid between t — t sc 
and t — 3t sc . As the pressure gradient is relieved by 
this expansion, the expansion slows and stops. The 
material remaining on the grid briefly falls back to- 
ward the center of the cluster before reaching hydro- 
static equilibrium. Because the computational box is 
not much larger than the cluster itself, at late times 
the density distribution in the 'vacuum' regions is not 
quite spherical. We use a larger box in the collision 
simulations to avoid this effect. 

The time required for the cluster to reach virial 
equilibrium, as determined using the method de- 
scribed in Section 5, was approximately (2-3)i sc in 
each test run. In order to control for this relaxation 
in the collision runs we place the cluster centers far 
enough apart initially that their cores do not interact 
until this amount of time has passed. 
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Fig. 3. — Angle-averaged density profiles in single- 
cluster test runs at t = 20t sc . The dotted line in- 
dicates the initial truncated Hubble-law profile, the 
dashed line indicates the final profile for run Tl, the 
dot-dashed line indicates run T2, and the solid line 
indicates run T3. 

Figure 3 shows angle-averaged density profiles for 



the three test runs at t = 20i sc , well after each has 
reached a steady state. We calculate each average 
profile by linearly interpolating from the computa- 
tional grid along 1000 randomly chosen lines emanat- 
ing from the center of mass; each line is uniformly 
sampled in radius out to the closest external grid 
boundary. The percentage deviation from average of 
the minimum and maximum density as a function of 
radius is plotted for each run in Figure 4. In run Tl, 
for which the central zone size was somewhat larger 
than r c , the final profile differs substantially from the 
initial profile, even in the outer regions of the cluster 
(2r c < r < R), where the density is almost a fac- 
tor of two larger than in the highest-resolution run. 
The density profile outside the cluster has not even 
converged, although its logarithmic slope has. In run 
T2, which used a central zone size of about 0.6r c , 
the outer regions of the cluster have converged to the 
high- resolution result, but the density at r « 2r c is 
overestimated by a factor of 1.2, and the region inside 
one core radius is low by a factor of 1.7. Although 
this resolution does not allow study of scales smaller 
than r c , it appears to be needed in order for the outer 
regions of the cluster to be correct. In both run Tl 
and run T2 the effect of poor central resolution shows 
up as an overestimate of the density in the zones just 
outside the center and an underestimate in the central 
zones. This effect is most likely due to errors in the 
central pressure and potential gradients arising from 
the use of a Cartesian grid, as can be seen in Figure 
4. The density asymmetry in runs Tl and T2 reaches 
its maximum in the zones just outside the center; for 
run Tl this maximum is about 35%, and for run T2 it 
is about 12%. Outside the cluster the asymmetry in- 
creases to more than 50% in each of the runs because 
of the finite box size effects mentioned earlier. 

The density profile in run T3, which used a zone 
size equal to one-fifth of a core radius, agrees with the 
64 3 result outside r « 2r c , but increases smoothly and 
monotonically to its central value, which is 20-30% 
smaller than the initial central value. This decrease 
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Fig. 4. — Fractional deviation from the average of 
the minimum and maximum density as a function of 
radius in the single-cluster test runs at t = 20t sc . 

in the central density is to be expected because of the 
initial pressure-driven expansion of the cluster, the 
consequent loss of 15% of the initial mass through 
the grid boundary, and the lack of a radiative cool- 
ing mechanism. For 2r c < r < 8r c the density drops 
as r~ 31 , slightly more steeply than the initial r~ 3 
profile. Outside of the cluster the density drops as 
r~ 4,8 . The density asymmetry inside the cluster is 
everywhere less than 5%, and in the outer regions of 
the cluster it drops as low as 1%. Because of the aver- 
age density profile's agreement with the 64 3 result, its 
smoothness and monotonicity, and its spherical sym- 
metry, we regard the 128 3 result as being adequately 
converged. The central zone size used in our colli- 
sion runs is approximately one-third of a core radius, 
between the zone sizes in runs T2 and T3. 
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Fig. 5. — Angle- averaged temperature profiles in 
single-cluster test runs at t — 20t sc . 



In none of the test runs does the cluster remain 
isothermal. During the cluster's initial expansion, 
most of the adiabatic cooling occurs in its outer re- 
gions. Figure 5 shows the angle-averaged temperature 
profile in each run at t = 20t sc . The temperature at 
the edge of the cluster (r = R) has fallen by a factor 
of five in each run, and at the edge of the grid it is 
smaller again by a factor of two. The central temper- 
ature is close to its original value of unity in each case, 
but in runs Tl and T2 it is about 5-10% higher. This 
is consistent with the underestimated central density 
and resultant shallowness of the density profile near 
the center in these runs. In run Tl the temperature 
in the outer regions of the cluster is badly underesti- 
mated. As with the density profile, the temperature 
profile in run T3 increases smoothly and monotoni- 
cally to its central value, 0.98. Outside of r 3r c it 
agrees well with the result from run T2. 

We conclude from these observations that resolving 
the core region is of great importance for minimizing 
numerical effects due to the Cartesian geometry of 
our grid. These effects include an underestimate of 
the central density relative to the regions just outside 
the center, asymmetry in the density profile just out- 
side the center, and excessive heating of the central 
zones. The amount of material lost from the cluster 
during its expansion, and the slope of the density and 
temperature profiles outside the cluster, do not ap- 
pear to depend significantly on the core resolution, 
but poor core resolution does cause an overestimate 
of the density and a substantial underestimate of the 
temperature in the regions just outside the core. A 
zone size smaller than r c /2 for the region r < 3r c 
appears to be needed in order to avoid these effects. 

Although we regard this resolution as adequate for 
the study of our nonradiating cluster gas, when com- 
paring our results to real systems the bremsstrahlung 
cooling time must be kept in mind. For radii smaller 
than r c the initial radiative cooling timescale, 

^= 2 - 6 (i) (mj^Y 1 {i-f ' ^ (5) 

can be comparable to the sound crossing time. In or- 
der to correctly resolve scales smaller than this it is 
necessary to include radiative cooling, which in addi- 
tion to permitting the study of cooling flows would 
tend to stabilize the cluster against the spreading 
caused by errors in the pressure and potential gra- 
dients. 
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4. Results 



In this section we present results from three col- 
lision simulations using equal-mass clusters at differ- 
ent impact parameters: zero, five, and ten core radii. 
The three collisions were followed until t = 12t sc , 
well beyond the point at which the merger remnant 
reached equilibrium in each case. The different col- 
lisions all involve multiple interacting shocks whose 
behavior and effect on the overall progress of the col- 
lision change with increasing impact parameter. Here 
we consider this variation as revealed by snapshots of 
the density, temperature, and velocity fields, and by 
the projected X-ray surface brightness and tempera- 
ture fields, at several different timesteps. 

4.1. Head-on collision 

The first run we performed, run A, simulated a 
head-on collision (&i n it = 0) . In Figures 6a and 6b we 
present snapshots of the density, temperature, and 
velocity fields at six representative times during the 
simulation in the plane having a constant z-value of 
1/2 the box height. The first figure shows an 8.1 R x 
8.1R section of the complete plane, and the second 
shows the innermost 16% of this region. In both plots 
the abscissa is the x-axis, and the ordinate is the y- 
axis. 

These figures show several important features of 
the collision. Of particular interest is the contrast 
between the behavior of the cores and that of the 
outer regions of each cluster. In the first snapshot of 
each figure, corresponding to t = 2t sc , the outer re- 
gions of the clusters have collided, forming a pair of 
slow-moving shock waves which propagate along the 
collision axis in both directions. Because of the clus- 
ters' gravitational potential, these shocks move more 
rapidly in the regions away from the collision axis, 
giving them a pinched appearance at the center of 
the grid. In the region between the shocks, the tem- 
perature is higher than its pre-shock value by about 
a factor of three. The approaching cluster cores com- 
press this shock-heated gas, causing it to be expelled 
from the collision axis at about one-half the clusters' 
initial free-fall velocity. The speed of this material 
increases by a factor of two as it accelerates down the 
pressure gradient. 

As the cluster cores collide at the center of the 
grid just before t = 3t sc , they produce a region of 
greatly enhanced pressure which is out of equilibrium 
with the local gravitational potential. This results in 




-2.o log (r/r lnl[ ) 2.0 

Fig. 6a. — Density and temperature snapshots for 
run A in the xy-plane passing through the cluster 
centers. Units arc as described in the text, with 
R = M = Tjnit = t sc = 1. Density contours are 
separated by a factor of three; the dashed contour in- 
dicates logp = — 1. Shading indicates the logarithm 
of the temperature. Velocity arrows are drawn for ev- 
ery eighth cell, with the fiducial arrow at the corner 
of each plot representing v = 1. 
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Fig. 6b. — Detail of the innermost 16% of Figure 6a 
(the xy-plane) . Log density contours are spaced every 
0.25. Velocity arrows are drawn for every sixth cell; 
the fiducial arrows indicate v — 0.5. 



a rapidly expanding, roughly ellipsoidal shock wave 
(Mach ~ 5.5) which heats the surrounding gas, car- 
rying away from the core the excess thermal energy 
resulting from the collision. This shock, which was 
seen also in the simulations of Schindler and Miiller 
(1993) (who described it as 'lens-shaped'), propagates 
more slowly in the region which was preheated by the 
initial pair of shocks, leading to a wishbone-shaped 
distortion in the ellipsoidal shock front. This is most 
clearly visible in the plots in the second row of Figures 
6a and 6b, which depict the simulation at t = 4t sc and 
t = 5t sc . At this point it has slowed by about a factor 
of two from its initial speed relative to the unshockcd 
gas, but because of the temperature gradient its Mach 
number has increased to about 18. 

In front of the ellipsoidal shock wave the initial 
pair of roughly planar shocks has continued to move 
parallel to the collision axis. Behind the ellipsoidal 
shock front, however, these shocks have mostly been 
disrupted, leaving behind a disk of relatively cold gas 
perpendicular to the collision axis. By t — 5t sc this 
gas has begun to fall back onto the merger remnant. 
The gas outside the region preheated by the initial 
pair of shocks (but now inside the ellipsoidal shock) 
has been shocked into radial motion away from the 
center of the grid, but is decelerating due the gravi- 
tational pull of the (now combined) cluster cores. 

By t — 6t sc the ellipsoidal shock has left the grid, 
taking with it approximately 15% of the initial mass 
of the clusters. The material left in its wake has begun 
to fall back onto the merger remnant, creating a weak 
cylinder-shaped accretion shock. The density inside 
the accretion shock is settling into a roughly spherical 
distribution, with small fluctuations in the velocity 
field. The temperature at the center is slightly higher 
than in the regions just inside the accretion shock. 
As time passes this temperature difference decreases. 
The central density and temperature after t = 7t sc are 
roughly constant at about twice their initial values. 

In Figure 7 we present projected X-ray surface 
brightness contour maps (integrated over the ROSAT 
energy band, 0.1-2.4 keV) for the merging system 
at four of the representative times of Figure 6a, as 
viewed along the x-, y-, and z-axis. For the purpose 
of generating these maps we have assumed a cluster 
mass M = 4.4 x 10 14 Af o , radius R = 1.58 Mpc, sound 
crossing time t sc = 1.5 Gyr, and initial temperature 
Tinit = 4.5 x 10 7 K. The bremsstrahlung emissivity eg 
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is given by (Rybicki & Lightman 1979) 



where, for our choices of G and /x/fe, 
We define 

r n/T 

I{T) = / x-° A e- x dx , (8) 

JT a /T 

where T a = 0.1 kcV//i and T 6 = 2.4 kcV/h define the 
observing energy range. (With the ROSAT energy 
band, I(T) is roughly constant and equal to unity for 
5 x 10 6 K < T < 10 8 K.) As Schindler and Miiller 
noted for their simulations, here we see that a favor- 
able viewing angle is needed to detect the ongoing 
collision; we note also that one must observe the sys- 
tem within one sound-crossing time or so of the core 
collision to see any significant distortions in the X- 
ray isophotes, even for the most favorable viewing an- 
gles (along the y and z axes). During this period the 
isophotes in the center of the map have a distinctly 
oval appearance. However, despite the strength of 
the ellipsoidal shock, the density enhancement just 
behind it is not sufficient to make it even as luminous 
as 10~ 5 times the central luminosity. The accretion 
shock is also not visible. In the X-ray projections at 
t = 5t sc the cluster appears spherically symmetric and 
relaxed from all directions. 

Projected temperature maps, weighted by the X- 
ray emission in the ASCA band (1.5-11 keV), are 
shown as grayscales in Figure 7. To generate these 
maps we assumed the same values for M, R, t sc , 
and Xlnit as we used in generating the X-ray surface 
brightness maps. The temperature maps make it clear 
that, in order to see the shocks in this problem (ex- 
cept for the initial pair), one must look not at the 
X-ray surface brightness but at the X-ray tempera- 
ture; all of the shocks described earlier are visible in 
these maps. This is true partly for the following rea- 
son. The bremsstrahlung emissivity (6) is very sen- 
sitive to the density and only weakly sensitive to the 
temperature. However, the Rankinc-Hugoniot con- 
ditions, which relate the jumps in density, pressure, 
and velocity across a shock to the shock's Mach num- 
ber M relative to the pre-shock gas, limit the density 
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Fig. 7. — Projected maps for run A of ROSAT X-ray 
surface brightness (S x ) and ASCA emission-weighted 
temperature (T x ), viewed along the x, y, and z axes. 
The innermost 1.6-R of the grid is shown. S x (in units 
of Mt~ 3 sr _1 ) is represented by contours spaced by a 
factor of three, with the outermost contour indicating 
log S x — —3.5. T x is represented by shading. 
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jump across the shock to (Landau & Lifschitz 1987) 
(7 + 1)M 2 7 + 1 



P2 
Pi 



(7+ 1) + (7 - 1)(M 2 - 1) M^oo 7-I 



(9) 



Here pi is the pre-shock density, and P2 is the post- 
shock density. On the other hand the temperature 
jump can increase without bound: 



[( 7 +l) + 2 7 (M 2 -l)] 
( 7 +l) 2 M 2 X 

[( 7 + l) + ( 7 -l)(M 2 -l)] 

27(7-1) 



(10) 



Af^oo (7 + If 



-M 



Naively substituting the limiting density (9) and tem- 
perature (10) ratios into the expression (6) for the 
bremsstrahlung emissivity, we find that the limiting 
cmissivity enhancement across a strong shock is 



£ff.2 



(7+1) 



2.2 



£ff,i M^oo (27)°- 1 (7 - 1) 



-0.1 



(11) 



where 



J(Ti) 



I rT b /T t 
I{Tl) J T a /Ti 



JlO.6 _ rpO.6 

06rp7(7\) 



(12) 

is a weak function of the pre-shock temperature. 
Hence X-ray emissivity (even without projection ef- 
fects) is a much poorer detector of strong planar 
shocks than gas temperature. It should be noted, 
however, that this result is only valid at large dis- 
tances from the center of the cluster potential. 

4.2. 6i n it = 5r c 

Run B was the first off-center collision, with impact 
parameter 6; n it = 5r c . Figures 8a and 8b show snap- 
shots of the density, temperature, and velocity in this 
run for 100% and the innermost 16%, respectively, of 
the midplane perpendicular to the z-axis (the colli- 
sion plane). Once again the abscissa in these plots is 
the x-axis, and the ordinate is the j/-axis. Figure 8c 
gives a zoomed view of the midplane perpendicular to 
the y-axis; here the ordinate is the z-axis. 

In this run, as in the head-on case, we see a se- 
quence of initial planar shocks, followed by a fast el- 
lipsoidal shock as the cores collide, followed by an ac- 
cretion shock as material falls back onto the merger 
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Fig. 8a. — Slices of log density and log temperature 
for run B, taken perpendicular to the z-axis through 
the center of the box. Units, ranges, and contour 
spacing are as in Figure 6a. 
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remnant. However, the behavior of these shocks is 
more complex, and the presence of angular momen- 
tum in the system alters their roles somewhat. The 
cluster cores survive intact longer, and the merger 
remnant requires more time to reach equilibrium. 

The planar shocks produced by the clusters' initial 
interaction (seen at t — 2t sc ) are now oblique, mak- 
ing a ~ 72° angle with the collision axis. The tem- 
perature jump and outflow velocity in the post-shock 
region are comparable to the values in run A. 

As previously noted, in the off-center runs the clus- 
ter cores survive intact beyond the time at which they 
would have merged in the head-on case. The size of 
each core is determined roughly by the ram pressure 
of the gas through which it passes: the cluster re- 
mains intact where its own gas pressure exceeds the 
ram pressure. As we increase the impact parameter of 
the collision, the density p' of the gas through which 
each core must pass decreases rapidly, while the rela- 
tive speed v remains roughly the same. The core's gas 
pressure p decreases with radius in the same way that 
p' decreases with &i n j t . Therefore the radius at which 
p drops below the ram pressure (~ p'v 2 ) should be 
roughly proportional to 6; n it- The sizes of the cluster 
cores seen in runs B and C are consistent with this 
explanation. 

Because the planar shocks are now oblique they 
play a new dynamical role not found in the head- 
on collision. The cluster cores, defined as described 
above, carry most of the angular momentum and sur- 
vive past the point at which they would have coa- 
lesced in the head-on case. The cores are not dis- 
rupted by the planar shocks; instead each passes 
through the shock which propagated toward it, then 
drives the other shock ahead of it, so that the cen- 
tral part of each planar shock is twisted into a spiral 
pattern. As each core moves through the outer re- 
gions of the other cluster, it creates a strong, curved 
bow shock which propagates out of the center of the 
grid. Together these form a fast ellipsoidal shock 
front similar to (but weaker than) that created by 
the core collision in the head-on case. This shock 
serves to dissipate the excess kinetic energy of the 
cluster cores, enabling them to fall into orbit about 
one another. However, not all of the excess is carried 
away by this shock; the remainder is dissipated more 
slowly by what is left of the slower, originally planar 
shocks (t = 3t sc ). By turning the material passing 
through them, these spiral shocks gradually transfer 
the orbital angular momentum of each core to the 
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Fig. 8b. — Detail of the innermost 16% of Figure 8 
(the xy-plane). Units, ranges, and contour spacin 
are as in Figure 6b. 
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Fig. 8c. — Slices of log density and log temperature 
for run B, taken perpendicular to the y-axis through 
the center of the box. Detail of the innermost 16% of 
the grid is shown. Units, ranges, and contour spacing 
are as in Figure 6b. 



surrounding low-density gas. Shortly after t = At sc 
the cores lose enough angular momentum in this way 
to fall almost directly in toward each other, having 
completed nearly a full revolution. 

By t — 5t sc the ellipsoidal shock has nearly left 
the grid, and the gas in its wake has begun to fall 
back onto the merger remnant, forming an accretion 
shock. By t — 7t sc this accretion shock is fairly well- 
developed and lies roughly 1.5R from the center of 
the remnant. The remnant itself is rotating, and weak 
spiral shocks (the remains of the twisted planar shocks 
described earlier) can be seen; these are completing 
the redistribution of angular momentum out of the 
cluster cores. By t = 7t sc the merger remnant has 
relaxed to a roughly ellipsoidal shape and is rotating 
about the z-axis. A small amount of density sub- 
structure is still present in the innermost part of the 
remnant (r < r c ) at this time, although on this scale 
the core of the remnant appears to be isothermal. 

The projected ROSAT X-ray surface brightness for 
run B (contours in Figure 9) shows a clearly sepa- 
rated double-peak structure while the cluster cores 
are orbiting one another. The bending of the pla- 
nar shocks is also visible. After the cores collide 
at t ~ 4i sc , the bimodal structure disappears, but 
the X-ray isophotes continue to appear significantly 
distorted, forming a barlike structure which rotates 
about the z-axis. When viewed along the bar, the 
cluster appears roughly spherical, and when viewed 
off-axis in the plane of the bar's rotation, the clus- 
ter appears elliptical, with the major axes of the 
isophotes roughly parallel. It is only when the clus- 
ter is viewed from well outside the plane of the bar's 
rotation that the differing ellipticities and angles of 
orientation of the isophotes reveal the ongoing equili- 
bration of the cluster. However, this cluster continues 
to appear aspherical from some directions long after 
the cluster formed in the head-on case has ceased to 
show any detectable evidence of a collision. 

The projected ASCA-weighted temperature for this 
run (grayscale in Figure 9) shows much clearer ev- 
idence for an off-center collision, but after the cores 
have collided this also requires a favorable viewing an- 
gle. As in the head-on case, a hot bar-shaped struc- 
ture is present in the temperature maps. However, 
here it is perpendicular to the line connecting the 
cluster cores, not to the initial collision axis. 
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Fig. 9. — Projected maps for run B of ROSAT X-ray 
surface brightness (S x ) and ASCA emission-weighted 
temperature (T x ), viewed along the x, y, and z axes. 
The innermost 1.6R of the grid is shown. S x (in units 
of Mt~ 3 sr _1 ) is represented by contours spaced by a 
factor of three, with the outermost contour indicating 
log S x = —3.5. T x is represented by shading. 

4.3. 6 init = 10r c 

In Run C we used an impact parameter b- m it = 
10r c . For this run, Figures 10a and 10b show full and 
zoomed snapshots of the midplane perpendicular to 
the z-axis (the collision plane). Figure 10c gives a 
zoomed view of the midplane perpendicular to the y- 
axis. This run displays many of the same qualitative 
features as the 6; n i t = 5r c case. However, while the 
collision itself is more violent, and the merger remnant 
requires much longer to relax to a new equilibrium 
configuration, the release of kinetic energy is much 
more gradual. 

The initial interaction of the clusters' outer regions 
again forms a pair of planar shocks, now at a ~ 54° 
angle with the collision axis. Because of their greater 
obliqueness they do not begin to twist until much later 
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Fig. 10a. — Slices of log density and log temperature 
for run C, taken perpendicular to the z-axis through 
the center of the box. Units, ranges, and contour 
spacing are as in Figure 6a. 

than in run B, and they brake the cluster cores and 
transfer their angular momentum to the surrounding 
gas less efficiently. The ellipsoidal shock which allows 
the cores to remain bound to each other does not 
form until just after t = 3t sc . In this run it dissipates 
a smaller fraction of the cores' initial energy than in 
the other two cases. A larger fraction of this energy 
is dissipated by the spiral shocks which are driven by 
the cores as they orbit each other. The cores com- 
plete about 11/2 orbits in the process of spiralling in 
toward each other; they do not merge until t = 6.5i sc . 
The size of the cores themselves is much larger (diam- 
eter ~ 8r c vs. ~ 5r c for 6; n i t = 5r c ), supporting the 
idea that their size is determined by ram pressure. 
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Fig. 10b. — Detail of the innermost 16% of Figure 10a 
(the iry-plane). Units, ranges, and contour spacing 
are as in Figure 6b. 



Fig. 10c. — Slices of log density and log temperature 
for run C, taken perpendicular to the y-axis through 
the center of the box. Detail of the innermost 16% of 
the grid is shown. Units, ranges, and contour spacing 
are as in Figure 6b. 
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The projected ROSAT X-ray emission for this run 
(contours in Figure 11) shows significantly more ev- 
idence for an ongoing merger than in the other two 
runs. Bimodal structure is present in the X-ray maps 
well past t — 5i sc . In the projected ASCA emission- 
weighted temperature maps (grayscale in Figure 11), 
bimodal structure also persists past t — 5t sc . The bow 
shocks are easily visible at t = At sc and 5t sc . From 
above the plane of the bar's rotation, some evidence 
of spiral structure can be seen at 5t sc . 

5. Discussion 

In this section we compare the three collision runs 
described in the previous section in terms of the time 
required for each to relax to its final equilibrium state, 
the increase in X-ray luminosity during the collision, 
and the structure of the merger remnant. 

5.1. Equilibration time 

In hierarchical models of large-scale structure for- 
mation, clusters of galaxies form through a series of 
mergers followed by periods of relaxation. How re- 
laxed clusters look in such a model at any given time 
depends on the ratio of the equilibration time to the 
average interval between mergers. If this ratio is much 
larger than one, so that, on average, clusters have not 
completely assimilated the effects of one merger be- 
fore the next takes place, clusters will not appear to 
be in equilibrium. If the opposite is true, so that the 
equilibration time is much smaller than the merger 
interval, then clusters will usually appear to be viri- 
alized, even if the merger rate is quite high. The 
fraction of clusters which appear to be far from equi- 
librium therefore is telling us something about this 
ratio. 

The rate at which mergers take place depends on 
the assumed underlying cosmological model, while the 
average equilibration time consists of two factors: the 
equilibration time for mergers having particular val- 
ues of the impact parameter, relative velocity, and 
relative mass, and the distribution of mergers with 
these parameters as a function of cosmological pa- 
rameters. The first factor is presumably only weakly 
dependent on cosmology, since the relaxing region has 
detached itself from the general cosmic expansion; the 
second depends on the tidal torques generated in dif- 
ferent models. If we are to reason backward from the 
fraction of clusters which appear unrelaxed and the 
rate at which clusters are observed to merge to de- 
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Fig. 11. — Projected maps for run C of ROSAT X-ray 
surface brightness (S x ) and ASCA emission-weighted 
temperature (T x ), viewed along the x, y, and z axes. 
The innermost 1.6-R of the grid is shown. S x (in units 
of Mi~ 3 sr _1 ) is represented by contours spaced by a 
factor of three, with the outermost contour indicating 
log S x — —3.5. T x is represented by shading. 
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ducc the underlying cosmological model (as specified 
by Clo, A, etc.), we need to understand the physical 
factors influencing the equilibration rate for collisions 
with given parameters. 

While our present simulations concern themselves 
only with the behavior of the gaseous component of 
galaxy clusters, thus ignoring the effect on the equi- 
libration rate of energy transfer between the gas and 
the dark matter, we can nevertheless draw some con- 
clusions regarding the rate of equilibration after an 
off-center merger. 

Physically motivated definitions for the equilibra- 
tion time will differ markedly from those based on 
observational criteria; furthermore, physically moti- 
vated definitions based on average quantities will dif- 
fer from those based on the actual presence of sub- 
structure (regardless of whether this is observable or 
not). In general we find that observational crite- 
ria (such as the projected luminosity and tempera- 
ture maps described in the previous section) lead to 
shorter equilibration times than criteria based on the 
degree of virialization of the merger remnant. 

From an observational standpoint the time re- 
quired in our simulations for the merging clusters to 
reach a new equilibrium state is short. The projected 
X-ray surface brightness (Figures 7, 9, and 11) shows 
little evidence of the shocks which are produced in 
each run. The high density of material in the cluster 
cores' gravitational wells accounts for most of the X- 
ray emission, and what evidence exists in the X-ray 
maps for an ongoing merger, such as bimodality, is 
only present when the cores are well-separated in the 
plane of the sky, or when they have not yet fully co- 
alesced. This is partly due to projection effects and 
partly because the bremsstrahlung emissivity poorly 
delineates planar shocks. 

In the head-on case bimodality is only present be- 
fore the collision; after the collision the anisotropy 
caused by material falling back onto the remnant per- 
pendicular to the collision axis is only visible for about 
one sound crossing time. In the case of 6; n it = 5r c 
the cores do not coalesce as rapidly; consequently 
they are distinguishable both before the collision and 
for more than one sound crossing time after the el- 
lipsoidal shock front forms. During this time they 
are orbiting one another, and they only appear well- 
separated if viewed from above the orbital plane or if, 
by chance, they are viewed at a point in their orbit 
when neither lies in front of the other. However, from 
all directions the X-ray isophotes remain fairly ellipti- 



cal, with an axial ratio of approximately 2:1, until sev- 
eral sound crossing times have passed. When viewed 
in the collision plane the long axes of the isophotes 
are all oriented in the same direction, parallel to the 
collision plane. When viewed from above the colli- 
sion plane, their axes are not aligned; the innermost 
contours are advanced in the direction of rotation rel- 
ative to the outermost ones. Initially the difference in 
angle is as large as 90°, but this decreases with time. 
By t — 5t sc the difference has decreased to 45°, and at 
late times the remnant appears spherically symmet- 
ric. Thus the rotation of the merger remnant in run B 
does not appear to significantly affect the shape of the 
X-ray isophotes at late times. The case of 6; n it = 10r c 
displays many of the same characteristics as that of 
&init = 5r c prior to the final merging of the cluster 
cores, except that the cores' inspiral requires more 
than three crossing times, allowing bimodal structure 
to be visible substantially longer than in runs A and 
B. 

The virial theorem provides a useful physical cri- 
terion for determining when the system is close to 
equilibrium. Whether the system is relaxed or not, 
the virial theorem requires that 



2 dt 



2 = 2T + W + 3(7 -1)E/- (13) 



(Chandrasekhar 1961), where I is the moment of in- 
ertia and T, W, and U are the total kinetic, potential, 
and thermal energies, respectively: 

I = \j d 3 ^(x)|x| 2 
T = \f d 3 xp(x)|v(x)| 2 

W/'^w (14) 

d 3 xp(x) . 
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The volume integrals in these definitions arc taken 
over the entire grid. When the system has reached a 
steady state, the left-hand side of equation (13) will be 
zero. We can approximate this criterion for practical 
purposes by dividing equation (13) by W: 



1 d 2 I 




2W dt 2 





„ T -, 3(7 -1W 

2 h 1 + — — 

W W 



< e < 1 . (15) 



When this criterion is satisfied for some suitably cho- 
sen value of e, the system will be close to equilibrium. 
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In Figure 12 we have plotted the second expression in 
equation (15) as a function of time for the three runs. 
We set e to 0.02, roughly the level of round-off error 
in our determination of T, U , and W, and roughly 
the point at which the virial parameter becomes con- 
stant in these runs. Although during the period fol- 
lowing the collision the behavior of the energy in run 
A differs substantially from that in run B, both runs 
reach equilibrium according to our criterion approx- 
imately 5t sc after the initial core interaction (repre- 
sented by the initial peak for each curve). Run C, in 
contrast, requires more than 6t sc to reach equilibrium 
after the initial core interaction. The virial parameter 
in this case passes through several maxima, reflecting 
the strength of the second and subsequent core inter- 
actions. 
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Fig. 12. — Factor by which the colliding clusters de- 
part from virial equilibrium as a function of time for 
the three collision runs. 

Using observational criteria to determine equilibra- 
tion time, we see that off-center collisions can take 
substantially longer than head-on collisions to reach 
a new equilibrium state. The virial equilibration time 
is less sensitive to the impact parameter but is at least 
twice the observational equilibration time. Given typ- 
ical sound crossing times t sc ~ 1 Gyr, both the obser- 
vational and the virial equilibration times can be a 
substantial fraction of the age of the universe. 

Adding collisionless dark matter to our simulations 
would tend to increase the equilibration time, since 
the dark matter does not experience shocks and so 
must dissipate its initial kinetic energy through vi- 
olent relaxation and drag due to the baryonic con- 
tribution to the potential. Adding radiative cooling 
would tend to decrease the equilibration time some- 



what by permitting the gas in the cluster cores to ra- 
diate away some of the kinetic energy in excess of the 
potential energy available to bind the cores, leading 
to weaker shocks during the collision. However, the 
cooling time is typically very long outside the cluster 
cores, so cooling should not affect the dynamics of 
these regions very much. In the central region where 
the cooling time is relatively short, the timescale for 
production of the ellipsoidal shock, which dissipates 
most of the cores' energy, is also short. Therefore it is 
likely that radiative cooling will have a much smaller 
effect on the equilibration timescale than it will have, 
for instance, on the luminosity of the clusters during 
the collision or the density and temperature profiles 
of the merger remnant. 

5.2. Brightening during the collision 

During the course of each collision the X-ray lu- 
minosity of the system varies significantly. The effect 
of changing the impact parameter on this luminosity 
evolution can be seen in Figure 13, which shows as a 
function of time the total energy loss rate, integrated 
over all frequencies, due to bremsstrahlung emission 
in each of the collision runs. This is just the total 
X-ray luminosity as a function of time: 

ijk 

(16) 

Before t = 1M SC the luminosity drops by about one- 
third due to the adiabatic cooling related to our initial 
conditions, as discussed in Section 3. Near t = 3t sc , 
when the ellipsoidal shock forms in each run, the ra- 
diative energy loss rate increases briefly by as much as 
a factor of 50, then drops to a new constant level 4- 
5 times its lowest pre-collision value. The height and 
width of the luminosity peak indicate the amount and 
duration of the compression experienced by the intra- 
cluster gas as the cluster cores approach one another. 
Thus in the head-on case we see one tall peak lasting 
about t sc /2, while in run B the initial peak is about 
one-third as tall and is followed by several smaller 
peaks as the cluster cores spiral inward toward one 
another. In run C the cores orbit several times before 
coalescing; hence we see several peaks, with the lu- 
minosity following each peak gradually increasing to 
roughly the same final level as in runs A and B. Be- 
cause of the long duration of the first orbit in run C, 
the second and subsequent luminosity peaks in this 
case appear about 3t sc after the first peak. 
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Fig. 13. — Radiative energy loss rates due to 
bremsstrahlung emission for the three collision runs 
(were cooling included). 



Because we do not include radiative cooling in 
our simulations, our luminosity enhancements are 
much greater than, for example, in the simulations by 
Schindlcr and Muller (1993). Nevertheless the results 
in Figure 13 provide information about the strength of 
the shocks produced in each run which should carry 
over to simulations which do include cooling. Our 
observations suggest that magnitude-limited surveys 
of galaxy clusters will be subject to a selection ef- 
fect based on merging activity, and in particular to 
one based on impact parameter. Head-on mergers 
produce stronger shocks and are more luminous than 
off-center mergers. Mergers are much more visible 
for a relatively short time while the cluster cores are 
first interacting than either before or after this time. 
Because mergers in progress can be much more lumi- 
nous than quiescent clusters, we may be overestimat- 
ing the amount of merging activity at high redshift. 
Alternatively, since projection effects can hide isopho- 
tal evidence for mergers at most orientations, we may 
observe merging systems as highly luminous, appar- 
ently relaxed clusters. 

5.3. Structure of the merger remnant 

For all three runs we were able to follow the 
progress of the collision until well after the merger 
remnant became virialized. In this section we com- 
pare the angle-averaged density and temperature pro- 
files in each run at t = 12£ sc . We also examine the 
rotation of the remnant in runs B and C at t — 12t sc . 

In Figure 14 we plot using solid lines the angle- 
averaged density profiles for runs A, B, and C at t = 



12t sc . To determine the effect of core heating on the 
structure of the merger remnant, we have fitted the 
average profiles for r < 2R using two models: the 
/3-model (Cavaliere & Fusco-Femiano 1976), 



Po 



[1 + (r/r c ) 
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(17) 



and the analytical approximation to the de Vau- 
couleurs (1948) model for elliptical galaxies described 
by Hernquist (1990), 
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(18) 



In Figure 14 the best-fit /3-modcl profiles are shown 
as dashed lines, and the best-fit Hernquist profiles 
are shown as dot-dashed lines. The parameters cor- 
responding to these fits appear in Table 2. Note that 
the initial density profile for each cluster (equation 1) 
corresponds to a /3-model profile with (3 = 1. 

Strictly speaking, the /3-model is intended to de- 
scribe the distribution of gas which is in equilibrium 
with a potential determined by a collisionless mat- 
ter component, such as galaxies or dark matter; the 
value of (3 then gives the ratio of the velocity disper- 
sion of the collisionless component to the temperature 
of the gas. Here we have no collisionless component, 
so (3 serves more as an indicator of deviations from 
an isothermal distribution and as a measure of dif- 
ferences in core radius and central density among the 
runs. In particular, we expect shock heating during 
a collision to lead to a merger remnant with a core 
radius larger than that of the original clusters. 

The Hernquist model, in contrast to the /3-model, 
diverges as r — > 0, and at large r it drops as r~ 4 , 
whereas fitting with the /3-model allows the slope at 
large r to vary. Because of our grid resolution and 
initial cluster profiles we do not expect the Hernquist 
model to give a good fit near r = 0, but fitting with 
it does allow us to check the asymptotic density slope 
independently of the shape of the density profile at 
small r. The asymptotic slope is of interest because 
A-body simulations of halo formation in cold dark 
matter (CDM) universes (proceeding, as it does in 
such models, through mergers) have yielded differing 
values for the slope. Methods neglecting infall from 
large scales, such as that used by Dubinski and Carl- 
berg (1991) in studying halo formation in CDM and 
that used by Hernquist (1992, 1993) in simulations of 
galaxy mergers, produce remnants with asymptotic 
profiles proportional to r~ 4 . Others (e. g. Navarro, 
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Table 2: Best- fit /3-model and Hernquist model parameters for collision runs at t = 12t sc . The first line gives 

/3-model parameter values for each cluster at t = 0. Note that the best-fit values for run C do not provide a 
satisfactory fit. 

/3-model Hernquist model 

Run p /MR- 3 r c /R (3 M H /M a H /R 

original 30.5 0.111 1.00 — — 

A 48.2 0.142 1.27 1.60 0.184 

B 60.6 0.119 1.19 1.69 0.205 

C 56.9 0.0904 1.02 1.74 0.304 
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Fig. 14. — Angle-averaged density profiles, taken 
about the center of mass, for the merger remnants 
in runs A, B, and C at t — 12i sc . The average profiles 
are plotted as solid lines. Best-fit model profiles de- 
scribed in the text include the /3-model (dashed lines) 
and the Hernquist model (dot-dashed lines). 



Frenk, & White 1996) have argued for an r -3 depen- 
dence on the basis of simulations which do include 
infall. Our calculations differ from these in that we 
only include collisional matter and we impose an ini- 
tial profile with r~ 3 asymptotic behavior. However, 
like Dubinski and Carlberg we implement a vacuum 
boundary, neglecting infall. If hydrodynamical effects 
are weak at large radius, we may therefore be able to 
check their results. In our single-cluster tests (Section 
3) we find that the asymptotic slope resulting from 
our artificial density cutoff at r — R (which would 
not be needed if we could include infall) is r 8 , much 
steeper than r~ 3 or r~ 4 , whereas the slope inside the 
cutoff radius remains close to r~ 3 . We may there- 
fore expect that any change to an r~ 4 profile in our 
simulations is an effect of the collision and not of the 
initial conditions. 

In examining Figure 14 and Table 2, we find that 
in none of the runs is the merger remnant isothermal; 
instead, the density profile in each case drops off more 
steeply than the initial profile at large radii. In runs 
A and B the profile for r > OAR is well-fit by the 
/3-model with similar values of (3. For r < OAR the 
/3-model produces a slightly poorer fit. As expected, 
the Hernquist model gives a poor fit for r < 0.2R; 
but at large radii both the /3-model and the Hern- 
quist model give good fits, with the best-fit value of 
[3 consistent with asymptotic slopes of 4.0 and 3.8 
for runs A and B, respectively, in agreement with the 
results of Hernquist (1992, 1993) and Dubinski and 
Carlberg (1991). For run C we were unable to obtain 
a statistically reasonable fit with either function; the 
density profile in this case is characterized by a large 
core (p drops to 1/2 its central value near r = 0.16i?), 
and while there is some evidence of power-law behav- 
ior outside the core, the logarithmic slope appears to 
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converge slowly. In all cases the density profile is con- 
sistent with an increased core radius. In runs A and 
B we see a large increase of the central density of the 
merger remnant in comparison with the central densi- 
ties in the original clusters, while in run C the central 
density is slightly smaller than the initial value. 

These results suggest that collisions with modestly 
varying impact parameter will produce merger rem- 
nants with substantially similar density profiles vary- 
ing as r~ 4 (if infall is neglected) at large radii, while 
the varying amount of shock heating in such collisions 
will primarily manifest itself in modest differences in 
the core radius. Most of the mass difference between 
the final merger remnant and each of the initial clus- 
ters goes into a large increase in the central density. 
A very large impact parameter, on the other hand, 
results in a significantly larger core with a more com- 
plicated asymptotic profile. In such cases the transfer 
of angular momentum out of the remnant core by spi- 
ral shocks (as discussed in Section 4) may not be ef- 
ficient enough to produce a remnant whose structure 
is independent of 6; n it- 



Run A (6 ini[ =0) 
Run B (b lnl , = 5r c ) 
Run C (6 Ml =10r c ) 



r / R 

Fig. 15. — Angle-averaged temperature profiles, 
taken about the center of mass, for the merger rem- 
nants in runs A, B, and C at t = 12t sc . 

The angle-averaged temperature profiles for the 
three runs at t = 12t sc , plotted in Figure 15, are con- 
sistent with these conclusions. For r < 2R the tem- 
perature profiles in runs A and B are virtually iden- 
tical, with run C slightly hotter than the other two 
for r < 0.1/? and substantially cooler outside this ra- 
dius. Inside r = 0.2i? the temperature in runs A and 
B is roughly constant, dropping from a central value 
of 2T; n it to about 1.5T in i t . Outside this radius the 



temperature drops approximately as r . It is inter- 
esting that the merger remnants in runs A and B are 
so similar; in the latter case the gas inside r = 0.3i? is 
rotating about the center of mass at an azimuthal ve- 
locity of about O.Si?^" 1 , but rotation in this case does 
not contribute much to its support against gravity. 
For 6i n it = 10r c , however, the rotation of the remnant 
has a very significant effect on both the density and 
temperature profiles, suggesting that if such extreme 
collisions occur at all, some evidence for them will be 
present in the structure of the merger remnant. 

6. Conclusions 

We have used a new hydrodynamical code based 
on PPM and a multigrid isolated potential solver to 
study the behavior of intracluster gas in controlled 
collisions between equal-mass galaxy clusters at dif- 
ferent impact parameters. A summary of our find- 
ings follows. Because of the lack of a collisionless 
component and of radiative cooling in our calcula- 
tions, as well as the simplified initial conditions, these 
conclusions should not be applied to specific clusters 
or cosmological models. Rather the physical insight 
suggested by these calculations concerning the behav- 
ior of self-gravitating gas spheres should be used to 
aid in the interpretation of more detailed calculations 
which include dark matter and radiation or which use 
more realistic initial conditions. The isolated poten- 
tial solver we have developed for these calculations is 
an important new component of a code we will use to 
perform more realistic calculations in the near future. 

Three main sets of shocks affect the progress of off- 
center mergers between clusters of equal mass. A pair 
of roughly planar shocks forms through the compres- 
sion of the clusters' outer regions as they approach one 
another. These shocks bound a disk-shaped outflow 
at a varying angle to the collision axis. In off-center 
collisions, these planar shocks become twisted into 
a spiral pattern through interaction with the cluster 
cores. The ram pressure acting on each cluster deter- 
mines the size of the cluster's core. The cores drive a 
strong, roughly ellipsoidal shock into the outer regions 
of the clusters, dissipating enough energy to permit 
themselves to fall into orbit about one another. In the 
wake of this shock the cluster cores drive two curved 
bow shocks which gradually transfer the angular mo- 
mentum of the cores to the surrounding gas as they 
spiral inward. After the cores collide, these spiral 
shocks continue to redistribute the angular momen- 
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turn and energy of the cores within the merger rem- 
nant. Material from the outer regions, heated and 
driven outward by the ellipsoidal shock front, then 
begins to fall back onto the merger remnant, creating 
a weak accretion shock. The Mach numbers of all of 
the shocks are modest. While the ellipsoidal and ac- 
cretion shocks have been seen in previous calculations 
of head-on mergers between clusters of galaxies (e. g. 
Schindlcr & Muller 1993), we only observe the spiral 
shocks in off-center mergers. We do not observe a 
high- velocity, ordered flow like that seen in the calcu- 
lations of Roettiger et al. (1993). This is most likely 
because of the low initial velocity and the lack of dark 
matter in our calculations. 

We find that both a physical criterion for equi- 
librium based on the virial theorem and the pres- 
ence of readily observable bimodal or elliptical struc- 
ture in the projected X-ray surface brightness and 
emission-weighted temperature lead to equilibration 
times which are a significant fraction of the age of 
the universe, given typical sound-crossing times of 
~ 1 Gyr. The virial equilibration time can be as 
large as 5-6 crossing times, while isophotal structure 
is erased within 1-2 crossing times after the initial 
core interaction for &j n j t = and b m i t — 5r c . For 
frinit = 10?*c the X-ray surface brightness shows ev- 
idence of structure for several crossing times due to 
the protracted inspiral of the cluster cores in this case. 
The head-on collision produces a much larger increase 
in brightness during the collision than the other two 
cases, but it is more short-lived; the luminosity en- 
hancement due to the collision lasts about one-half of 
a crossing time. 

The angle-averaged density and temperature pro- 
files of the merger remnant show very little variation 
between a head-on collision and an off-center collision 
with an impact parameter equal to five times the core 
radius. The core radius increases slightly from its ini- 
tial value due to shock heating during the collision, 
and the central density and temperature finish with 
values about twice their initial values. The central 
0.2R of the remnant is roughly isothermal; outside 
this radius the temperature drops as 1/r. Outside the 
core the density drops steeply, with a (3 value of 1.2 
— 1.3, corresponding to an asymptotic behavior closer 
to r~ 4 than to r~ 3 . Although the merger remnant 
rotates in the off-center case, the remnant does not 
appear to be rotationally supported. For an impact 
parameter of ten core radii, however, the remnant is 
characterized by a large core, a low central density, 



and a complicated profile at large radius, suggesting 
that the outward transfer of angular momentum by 
spiral shocks is inefficient in this case. 
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